{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using Fermi-LAT All-sky Weekly Files\n",
    "\n",
    "Some types of data analyses lend themselves to the use of all-sky data files, rather than repeatedly downloading overlapping regions of sky for repeat analysis. Examples are:\n",
    "\n",
    "* Light curves and spectra for a large number of sources, whether galactic, extragalactic, or both\n",
    "* Characterization of large-scale diffuse emission\n",
    "* Analysis of moving sources, like the Sun\n",
    "* Searches for transient events (like GRBs) that have an isotropic distribution\n",
    "\n",
    "This analysis thread describes how to generate an all-sky data set for use in other types of analysis. We will then create an exposure-corrected all-sky image as an example of how to use the weekly files."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Download the weekly LAT files\n",
    "\n",
    "The weekly LAT data files are available from the FSSC [FTP Server](https://heasarc.gsfc.nasa.gov/FTP/fermi/data/lat/). These files contain the LAT data acquired during each week of Fermi's science mission. (Mission weeks run Thursday through Wednesday, UTC time.) The file for the current mission week will also be listed, even though the week may not yet be complete.\n",
    "\n",
    "The current mission week file is regenerated by the FSSC after every data delivery from the LAT instrument team. This guarantees that the FTP Server and the LAT data server contain the same data, and both are up-to-date.\n",
    "\n",
    "You can download the LAT weekly files one at a time from the FTP site, or you can use `wget` to download a full set of files.\n",
    "\n",
    "Once you have the current set of files, you can use the same `wget` command (from the same directory) to download only the newest files. In the event that the LAT team reprocesses data for the full mission, the filenames will change, and `wget` will download a full set of the new, reprocessed, data.\n",
    "\n",
    ">**Note**: This may take a while."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget -m -P . -nH --cut-dirs=4 -np -e robots=off https://heasarc.gsfc.nasa.gov/FTP/fermi/data/lat/weekly/photon/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "LAT data has been classified based on the quality of the event reconstruction. The classification can then be used to reduce background and improve performance by filtering out poorly reconstructed events.\n",
    "\n",
    "With the advent of Pass 8 data reconstruction, there are a myriad of various event classes that are tailored for specific types of analysis. The user should see the [Cicerone](https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Data_Exploration/Data_preparation.html) pages for a full description of the various event classes (and types) available in Pass 8.\n",
    "\n",
    "For most analyses, the user will be interested in `SOURCE` class photons (`evclass=128`). However, analyses that are intrinsically signal limited require very low background. For example, a user studying large scale, diffuse structures may be interested in using the `CLEAN` (`evclass=256`) or even `ULTRACLEAN` (`evclass=512`) classes.\n",
    "\n",
    "There are two types of weekly files available; photon files and extended files. The photon weekly files contain all the SOURCE class and are usable for most analyses. The extended files contain the same events as the photon files, but also contain additional reconstruction information about each photon that may be useful if you are trying to characterize the quality of a particular signal, such as for dark matter line searches.\n",
    "\n",
    ">**NOTE**: For analysis of short timescale events (for example, GRBs) where the flux of the event is much brighter than the accrued background, the user may be interested in using the `TRANSIENT010` (`evclass=64`) event class. If you are interested in investigating transient events within LAT data, you cannot use the weekly files, as they have been prefiltered to contain SOURCE class photons. There are no weekly files that contain the `TRANSIENT` class events.\n",
    ">\n",
    ">In order to generate an all-sky dataset of `TRANSIENT` class events, you will need to perform a grid of [data server](https://fermi.gsfc.nasa.gov/cgi-bin/ssc/LAT/LATDataQuery.cgi) queries with the maximum radius of 60 degrees. Once you have downloaded data from the full sky, you will need to combine the files and filter out duplicate events. We recommend using the [FTOOLS](https://heasarc.gsfc.nasa.gov/lheasoft/ftools/ftools_menu.html) suite to perform these tasks.\n",
    "\n",
    "Run this command to download the weekly extended files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget -m -P . -nH --cut-dirs=4 -np -e robots=off https://heasarc.gsfc.nasa.gov/FTP/fermi/data/lat/weekly/extended/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to do analysis for the full mission, you will also need the mission-long spacecraft file. This file contains entries at 30-second intervals during all periods of active data-taking by the LAT instrument. As a result, this file can be rather large.\n",
    "\n",
    "There are also weekly spacecraft files that contain the same information. However, some users have found the pre-generated full mission file to be more robust during analysis than a concatenated set of weekly files.\n",
    "\n",
    "We recommend re-downloading the full mission file when you update your all-sky data set.\n",
    "\n",
    "Run this command to download the mission spacecraft file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget -m -P . -nH --cut-dirs=4 -np -e robots=off https://heasarc.gsfc.nasa.gov/FTP/fermi/data/lat/mission/spacecraft/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Combine the data files\n",
    "\n",
    "You should now have a full set of weekly files and a mission-long spacecraft file. The next step is to combine the weekly files into a single file. We will use the [gtselect](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtselect.txt) tool which can also be used for filtering data. But we want to combine the full data set without removing any events, so we will use slightly different parameters.\n",
    "\n",
    "To combine files, we first create a file list:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls ./weekly/photon/lat_photon_weekly* > ./data/filelist.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat ./data/filelist.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you look at filelist.txt, you will see it is simply a list of all the files in the directory that start with the string `lat_photon_weekly`. You will also want to clear any parameters in gtselect from previous analysis that could inadvertently remove events from the dataset.\n",
    "\n",
    "Clear the **gtselect** parameter file back to defaults with the command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!punlearn gtselect"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you're ready to run [gtselect](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtselect.txt) to combine the data files. Here we will use `INDEF` so that gtselect will not make a selection on event class or event type. See the **gtselect** [help file](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtselect.txt) for a description of the evclass and evtype cuts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtselect evclass=INDEF evtype=INDEF\n",
    "    @./data/filelist.txt\n",
    "    ./data/lat_alldata.fits\n",
    "    0\n",
    "    0\n",
    "    180\n",
    "    INDEF\n",
    "    INDEF\n",
    "    30\n",
    "    1000000\n",
    "    180"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combining the files may take a while.\n",
    "\n",
    "Now that you have mission-long data and spacecraft files, you are ready to generate AGN light curves, search for flaring sources, or perform any other type of analysis that requires using data from many regions on the sky. A simple example is described below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using the LAT data to produce an all-sky image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 1: Remove limb photons"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, filter the data for the proper event class and to remove the earth limb contamination.\n",
    "\n",
    "We will use the `SOURCE` class (`evclass=128`) as it has been tuned to balance statistics with background for long-duration point source analysis.\n",
    "\n",
    "Here we also include the `evtype=3` cut to use both front and back converting event types. We will remove contamination by limiting the reconstructed zenith angle to events at an angle of 90° or less. We will also improve the PSF by excluding events with reconstructed energies below 1 GeV.\n",
    "\n",
    "Run [gtselect](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtselect.txt) to filter the data and apply a zenith cut."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtselect\n",
    "    ./data/lat_alldata.fits\n",
    "    ./data/lat_source_zmax90_gt1gev.fits\n",
    "    0\n",
    "    0\n",
    "    180\n",
    "    INDEF\n",
    "    INDEF\n",
    "    1000\n",
    "    500000\n",
    "    90"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, correct the exposure for the events you filtered out.\n",
    "\n",
    "The way the Fermitools account for exposure is computed based on the Good Time Intervals (GTIs) recorded in the event file. Because we have eliminated some events, we need to update the GTIs for the filtered file. Even though we have used a zenith cut on the data, we will not correct for it here.\n",
    "\n",
    "For an all-sky analysis, an ROI-based zenith cut would eliminate the entire dataset. Instead, we will later correct the livetime calculation for the zenith angle cut made with [gtselect](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtselect.txt).\n",
    "\n",
    "The [gtmktime tool](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtmktime.txt) updates the GTI list based on events we eliminated in the previous step. It also provides the ability to filter on any of the parameters stored in the spacecraft pointing history file. We will use the most basic filter recommended by the LAT team."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtmktime\n",
    "    ./mission/spacecraft/lat_spacecraft_merged.fits\n",
    "    (DATA_QUAL>0)&&(LAT_CONFIG==1)\n",
    "    no\n",
    "    ./data/lat_source_zmax90_gt1gev.fits\n",
    "    ./data/lat_source_zmax90_gt1gev_gti.fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2: Bin the data\n",
    "\n",
    "Now you will bin the data in preparation for exposure correction. In order to correct the map for exposure, you have to account for the changing response with energy. To do this, bin the data into a counts cube using [gtbin](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtbin.txt), but with a single energy bin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtbin\n",
    "    CCUBE\n",
    "    ./data/lat_source_zmax90_gt1gev_gti.fits\n",
    "    ./data/lat_source_zmax90_gt1gev_ccube.fits\n",
    "    ./mission/spacecraft/lat_spacecraft_merged.fits\n",
    "    3600\n",
    "    1800\n",
    "    0.1\n",
    "    GAL\n",
    "    0\n",
    "    0\n",
    "    0\n",
    "    AIT\n",
    "    LOG\n",
    "    1000\n",
    "    500000\n",
    "    1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 3: Calculate exposure map\n",
    "\n",
    "Because Fermi data are sparse, the accumulated exposure time for any given point in the sky is based on the observatory's pointing history and the response of the instrument to events at various incidence angles and different energies.\n",
    "\n",
    "The LAT team has calculated how the instrument responds by modeling the entire instrument and then running a MonteCarlo simulation of many millions of events of different energies and directions. These simulations characterize how the gamma rays propagate through an ideal instrument, and the information is recorded in the Instrument Response Functions (IRF). Since the response depends on how much data filtering is performed, there are different IRFs for the different event classes.\n",
    "\n",
    "To see which IRFs are available from within the Fermitools, run the [gtirfs](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtirfs.txt) command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!gtirfs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The appropriate IRF for this data set is `P8R3_SOURCE_V2`.\n",
    "\n",
    "Calculating exposure for an all-sky map first requires calculating the instrument livetime for the entire sky (using [gtltcube](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtltcube.txt)) and then convolving the livetime with the IRF (using [gtexpmap2](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtexpmap.txt)).\n",
    "\n",
    "The pixel size for the livetime does not need to match the binned data, but you do need to provide the event file as [gtltcube](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtltcube.txt) uses the good time intervals to calculate the livetime at each position.\n",
    "\n",
    "You will also need to correct the livetime for the zenith angle cut made with [gtselect](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtselect.txt) earlier. To do this, use the `zmax` option on the command line and match the value used with [gtselect](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtselect.txt). This modifies the livetime calculation to account for the events you removed earlier.\n",
    "Here is an example of the input for [gtltcube](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtltcube.txt):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtltcube zmax=90\n",
    "    ./data/lat_source_zmax90_gt1gev_gti.fits\n",
    "    ./mission/spacecraft/lat_spacecraft_merged.fits\n",
    "    ./data/lat_source_zmax90_gt1gev_ltcube.fits\n",
    "    0.025\n",
    "    1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you can calculate the exposure at each point in the sky in each energy bin (we used only one bin) using [gtexpcube2](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtexpcube2.txt).\n",
    "\n",
    "You will need to enter the pixel geometry and energy binning to match what you used in the counts cube.\n",
    "\n",
    "Also, be sure to use `CENTER` for the energy layer reference. If you use `EDGE`, you will get an additional energy plane in the output file that causes complications later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtexpcube2\n",
    "    ./data/lat_source_zmax90_gt1gev_ltcube.fits\n",
    "    none\n",
    "    ./data/lat_source_zmax90_gt1gev_expcube1.fits\n",
    "    P8R3_SOURCE_V2\n",
    "    3600\n",
    "    1800\n",
    "    0.1\n",
    "    0\n",
    "    0\n",
    "    0\n",
    "    AIT\n",
    "    GAL\n",
    "    1000\n",
    "    500000\n",
    "    1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note**: You could use CALDB for the IRF, but if this doesn't work, try explicitly specifying the IRF. The appropriate IRF for this data set is `P8R3_SOURCE_V2`.\n",
    "\n",
    "If you open the exposure map with _ds9_, you will see that the image has two energy planes. The first plane is for 1-300 GeV, the second is for 300 GeV and above. We will only use the first image plane to correct our all-sky map."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 4: Correct the all-sky image for exposure and scale\n",
    "\n",
    "To complete the exposure correction, you must perform a set of arithmetic functions on the all-sky counts cube using the [FTOOLS](https://heasarc.nasa.gov/lheasoft/ftools/ftools_subpacks.html) **fextract**, **faith**, **fimgtrim**, and **fcarith**. First, extract the first of the two exposure planes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "fextract\n",
    "    ./data/lat_source_zmax90_gt1gev_expcube1.fits[0]\n",
    "    ./data/testfextract.fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, correct the value of each pixel for the exposure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "farith\n",
    "    ./data/lat_source_zmax90_gt1gev_ccube.fits\n",
    "    ./data/lat_source_zmax90_gt1gev_expcube1.fits\n",
    "    ./data/lat_source_zmax90_gt1gev_corrmap.fits\n",
    "    DIV"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then trim the pixels that were outside the Aitoff projection to zero so they don't affect your display."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "fimgtrim\n",
    "    ./data/lat_source_zmax90_gt1gev_corrmap.fits\n",
    "    0\n",
    "    0\n",
    "    INDEF\n",
    "    ./data/lat_source_zmax90_gt1gev_corrmap2.fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, scale the image so that the maximum pixel is equal to 255 (the typical maximum dynamic range for graphics).\n",
    "\n",
    "There are a number of ways to find the maximum pixel value. Here we have used _ds9_ to find a maximum pixel value of $4.143x10^{-8}$. This means our scale factor should be $255/4.143x10^{-8}$, or $6.155x10^9$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "fcarith\n",
    "    ./data/lat_source_zmax90_gt1gev_corrmap2.fits\n",
    "    6.155e9\n",
    "    ./data/lat_source_zmax90_gt1gev_corrmap3.fits\n",
    "    MUL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can look at the final product in _ds9_ using log scaling. This process is approximately the same as is used to create the LAT all-sky images released by the instrument team.\n",
    "\n",
    "This method can be used for different energy ranges, different event class selections, or even to compare different data sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
